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Taking advantage of a novel approach to the design of the orbital transfer optimization problem and 
advanced non-linear programming algorithms, several optimal transfer trajectories are found for problems 
with and without known analytic solutions. This method treats the fixed known gravitational constants as 
optimization variables in order to reduce the need for an advanced initial guess. Complex periodic orbits 
are targeted with very simple guesses and the ability to find optimal transfers in spite of these bad guesses 
is successfully demonstrated. Impulsive transfers are considered for orbits in both the 2-body frame as well 
as the circular restricted three-body problem (CRTBP). The results with this new approach demonstrate the 
potential for increasing robustness for all types of orbit transfer problems. 


Introduction 

Solving complex trajectory problems optimally 
is of vital importance to the success of future 
aerospace programs. For example, NASA is 
currently engaged in developing a new human 
spaceflight program to the Moon that places new 
complex constraints on the spacecraft trajectories 
there and back. This is only one example; other 
interplanetary missions are being developed and will 
lead to new programs that will continue to expand 
out into the universe. If we are to succeed in these 
endeavors, our ability to find optimal trajectories 
must expand with it. 

In the ever growing field of optimization 
research, algorithm robustness is the ultimate aim. 
One approach to achieve this robustness is to modify 
current optimization approaches to handle complex 
problems without requiring meticulous tuning to 
obtain a solution; in other words, construct an 
algorithm that does not require a feasible target 
solution as a starting point or advanced knowledge 
of the direction of the optimal solution. 

Through the ages, trajectory optimization 
research has been divided into two main categories. 
Both methodologies have been tested and 
extensively applied to specific orbital trajectory 
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problems. Indirect methods take advantage of 
calculus of variations techniques to solve two-point 
boundary problems. Direct methods transform the 
continuous problem into a parameter optimization 
problem which is solved numerically. Many 
approaches to the direct problem have been 
developed including implicit approximation of the 
states as well as explicit integration techniques, 
known as the direct shooting approach. 

The algorithm presented here employs the 
second type of direct trajectory optimization, but 
with a new twist. The equations of motions that 
dictate the behavior of bodies in space have fixed 
constants. These constants may be invariable, but 
that does not mean that they cannot be treated as 
optimization variables. The constants can be 
constrained to the fixed values, thus imparting a 
target error on known values and ensuring the final 
solution is feasible. 

p as an Optimization Parameter 

For orbit dynamicists, the primary constant that 
appears in almost all calculations is the gravitational 
parameter, determined individually for each massive 
body in the solar system. Imagine treating this 
parameter as a variable. In the simplest example, a 
two-body system, the force model is given by the 
following equation of motion: 
jtir 
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The dynamics in this system is controlled 
entirely by the value of this parameter. If the 
parameter is very large, like the sun, the force 
exerted on the spacecraft is equally large. But if the 
gravitational constant were very small this force is 



similarly smaller and less influential. In fact, if the 
parameter was small enough to be infinitesimally 
tiny, a force-free field is created. In this system, an 
optimal trajectory has an obvious structure: a 
straight line. 

It is this fact that the algorithm attempts to 
exploit. Start with an easily constructed straight line 
guess in a field free system and the algorithm adapts 
the solution to the unique quality of the problem by 
changing values of the gravitational parameter in the 
direction of the known constant. If the solution 
process from iteration to iteration is modified 
gradually to the desired value, the optimization 
routine may be able to handle very complex 
problems. 

In applying this technique to more complex 
systems, such as the circular restricted three-body 
problem (CRTBP), the gravitational parameter has a 
different meaning and a straight line guess is no 
longer precisely a solution at any phase of the 
problem. For example, in the CRTBP p is now a 
ratio of gravitational parameters, the secondary over 
the primary body. A value of zero simply eliminates 
the secondary body from the problem, the primary 
body is still present and thus a force-free field space 
can never exist. However, this fact does not 
necessarily discount the usefulness of treating the 
parameter as an optimization variable in the CRTBP 
frame. 

The CRTBP system is useful tool for mission 
design because of the presence of periodic orbits that 
subsist solely due to the interaction of the two 
bodies. If the gravitational parameter changes these 
periodic orbits breakdown. However, transfers to 
these orbits typically take advantage of the unstable 
manifolds associated with these orbits, so the 
breakdown that occurs with the changing parameter 
may be useful in zeroing in on an optimal solution. 

While a straight line guess is no longer valid it is 
still a useful starting point. In this instance, this is 
particularly made true if the number of independent 
segments is large. The segments can then be splayed 
out equally between the two targets and along a 
straight line. Propagating each segment does not 
deviate very significantly from this line initially. 

Direct Multiple Shooting Approach 

The segments in the previous section refer to the 
multiple shooting arcs in the formulation of the 
optimization problem. In this particular multiple 
shooting formulation every segment contains a set of 
8 independent variables including initial and final 
time and 3 components each of position and velocity 


at the initial time. Since the number of ideal 
segments might vary from problem to problem the 
algorithm was designed to handle any n number of 
them. Each segment thus has the following 
parameter structure: 

= [-To y t o z io vx i0 vy i0 vz, 0 t i0 t if \ 


From this parameter vector, the state is propagated 
to the final time. Continuity is then ensured not 
internally but by a constraint parameter vector. The 
state and time at the endpoints of the segments are 
ensured to be continuous by equality constraints and 
time from segment to segment is guaranteed to 
increase monotonically with inequality constraints: 
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Next, constant p needs to be included as a 
parameter. This variable is the same value for each 
segment, there is not a separate gravity constant for 
each segment: 


( 3 ) 
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Equations 2 and 3 together form the basic structure 
for the variables that are used in every formulation 
of the problem. 

The differences occur in the end-point target 
constraints. However, throughout this paper, the 
target constraints are the same. They are all position 
constraints (velocity differences at the endpoints are 
incorporated in the objective function). Assuming 


additional time like parameters ( r 0 and T f ) for the 


initial and final orbits, the constraints are given 
below: 
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Formal Problem Definition 

With these constraints and parameters 
established, it is possible to delineate the entire 
problem formulation. Since impulsive transfers are 
the focus of this algorithm, the initial and final Av’s 
constitute the two parts of the objective function: 

J = Av 0 + Av f (5) 

In terms of the parameter vector variables, this is 
given as follows: 

j = [(IT*,- - vx io ) 2 + (XV W - vy,o ) 2 + (yz M - VZ 10 ) 2 f ' 2 + 

r -i. (u) 

[( vx n f - vx ,a, ) 2 + (yy„f - yy ur ) 2 + ( vz „ f - vz , a , ) 2 J ' 

And finally the entire parameter vector and 
constraint vector combined: 



The length of the parameter vector is 8n + 3. 
The constraint vector has length 8n + 1 . 8n of these 
constraints are inequality while the rest are equality 
constraints. Since the number of total constraints 
does not exceed the length of the parameter vector 
this is a well-defined optimization problem. (The 
number of variables only has to exceed the number 
of equality and active inequality constraints, but 
noting the former difference guarantees this.) It is 
also worthwhile to note that the inequality 
constraints will only be active if the segment is 
reduced to be infinitesimally small (t i0 = t if ). All of 
the derivations in this paper are based solely on the 
parameter and constraint vectors above, but some of 


the results in the 2 body system use true anomalies 
in place of x’s for the final two parameters; the 
derivations of these gradients have not been 
provided here for brevity. 

Analytic Gradient Derivation 

In many instances, for the sake of expediency, it 
is logical to avoid the tedious derivation process 
necessary to obtain precise analytic derivatives and 
use a numerical finite difference method instead. 
This is not one of those instances. While numerical 
methods can provide very close approximations of 
the actual derivatives nothing can approach the 
accuracy of analytical gradients. And while for some 
simple problems or problems that have good initial 
guesses these derivatives are adequate, the goal of 
this formulation is to be robust, and thus analytic 
gradients are necessary to ensure a truly powerful 
algorithm. 

In addition, with the size of this problem 
potentially being very large depending on the 
number of segments selected, the importance of 
analytic gradients in terms of the efficiency and/or 
capability in solving a given problem becomes even 
more critical. With finite differences, each additional 
variable and constraint adds at least a couple of 
fimction evaluation steps (or several) to each 
iteration. With analytic derivatives, the relative time 
spent in an iteration step increases only as much as 
the derivative function increases in complexity. 

In general, most of the derivatives for this 
formulation are independent of the specific problem 
being solved. In this case, since the targets always 
contain position constraints, not even these end- 
point derivatives need to be re-derived. 

In order to obtain all of these gradients linear 
perturbation theory must be employed. This theory 
posits that if the state, delineated by the vector x , is 
perturbed along a path,T(/), at some instance t 
along that path, this variation can be related to a 
corresponding variation in the initial state at the 
initial time (or any other arbitrary time) by use of the 
following equation: 

fic(t) = ®(t,t 0 )fic(t 0 ) (8) 

Here the matrix &(t,to) refers to the state transition 
matrix (STM), which can be obtained by integrating 
the state propagation matrix F, along with the 
equations of motion. The state propagation matrix 
(SPM) is simply the partial derivative of the time 
derivative of the state equation vector with respect to 
the state vector. 



Augmented State Propagation Matrix 

In order to get the correct gradients for (i, the 
gravitational parameter must be added to the state 
vector A relation between a variation in p and the 
initial state of any particular segment is needed. 
Thus the augmented state vector becomes: 
frl 


The state propagation matrix is then: 
f dv dv dv 


f = ^ l = ^ = 


Everything in the above matrix, which depends 
directly on the force field propagated, has been 
derived before except for the last row and last 
column: the dependence on p. Since this algorithm is 
used to solve problems in both the 2-body and the 
CRTBP systems, the derivation for both will be 
given. In all derivations the time derivative rate of p 
is assumed to be zero. The derivation in the 2-body 
frame is very simple: 

djii _ djii _ djii _ dv _ ^ 
dr dv du du 

(11) 

dr r 

d/u r 3 

The augmented propagation matrix in the two- 
body frame is thus: 
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The same derivation in the CRTBP is a bit more 
complicated. The CRTBP force field equations of 
motion are as follows: 
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where: 

r n = \[x + fjj 2 +y 2 +z 2 }' 2 

r 23 = [(* - (i - ju)) 2 + y 2 + z 2 1 2 
As before: 

d/u _ dju _ dju _ dv _ 
dr dv dju dju 

Start by finding the partials of the r vectors: 
dr u _ x + /u 
d/u r u 
dr 23 _ (x + (!-//)) 
d/u r 23 

Then: 
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The entire propagation matrix for the CRTBP is thus: 
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Objective and Constraint Derivatives 

The usefulness of the state transition matrix 
found by propagating the SPM will be readily 
apparent after objective and constraint equations are 
differentiated with respect to the parameter vector. 

In this derivation the same partial derivative is 
repeated multiple times, so only the original partial 
will be shown here and assumed trivial to repeat. 
Start with c i7 (the first constraint ensures continuity, 
so the indices i and i -1 will be used for the 
constraints and state vectors respectively): 
dc d(x.„ —x„ n , ) dx„ ,,, 

UL ’' = v ,0 (H 3 LL = 2d22_ = _<r> (t t ) MO! 

dx 8x 8x 1,1 ('-D/ ’('-i> oJ 

^(1-1)0 l/JI '(i'-l)0 CM '(i-l)0 

The validity of using the state transition matrix 
goes back to linear perturbation theory. This can be 
applied in every instance that a partial of a state 
vector of a segment at a final time is needed with 
respect to an initial time. The next challenging 
gradient is the partial derivative with respect to the 
initial time: 

SCg _ d(*,-0 -*(,--!)/) _ foi-l)/ 

^(i-l)0 ^(f-l)0 ^(i-l)0 

The derivative at this point is not directly 
relatable to the STM, but the relationship between 
the differential and the time independent variation is 
known: 
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And can be substituted: 
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This can be simplified with the STM and another 
substitution of the differential (this time with x (i _i )0 )'. 
dc 

TT = i (^(i-i)/ ’ ^(i-i)o 1)0 (24) 

CT (i-l)0 

Now evaluate with respect to the final time: 

3c,. i 

TT A - = - vx Wf (25) 

Compared to the previous state, the partial of 
this constraint to the current state i is relatively 
trivial, it is simply the identity matrix. 

The importance of augmenting the SPM is 
evident in the partial with respect to p: 

dc a 4,0 -Viv) ... , 

dju~ dv ~ d M - 17((w)/ ’ ('-DO - 1 (26) 


The x’s only come into play with the end point 
constraints and are also trivial: 

d(*w~*i 0 ) = ( ) (27) 


S(x nf ~x tar ) 


= -VX tar {Tf) 


( 22 ) 



The full derivative matrix can now be created for n segments: 
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The potential for mistakes is obvious, so all of these gradients were verified by numerical approximations to 
within a desired tolerance. 


The objective gradients can be obtained in a 
similar manner. Recall that: 
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Then the objective derivatives are: 
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Objective derivatives (cont): 
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Primer Vector Theory 

The optimal control theory that is used to 
develop primer vector theory has been presented 
many times, so only a summary of the aspects used 
in this optimization scheme are given here. 

The theory that applies to this particular case is 
restricted to an objective function that sums 
impulsive velocity increments, in this case, at the 
endpoints. It has been shown, in both the two-body 
system and the CRTBP system that the primer 
vector and primer vector rate adhere to the same 



equations as the state variations. For the two-body 
system this means that: 
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Likewise for the CRTBP: 
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Where G c is not the same matrix as G. 

These equations make it possible to solve for the 
primer vector and the primer vector rate in terms of a 
given state and the STM that is used to transfer it: 
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This equation is useful if relationships for at least 
two out of the four quantities in the above equation 
are available. Fortunately, the primer vector is 
relatable to known values of the state by Lawden’s 
necessary conditions for an optimal trajectory: 
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With these quantities known, the above equation is 
now a two point boundary value problem, which can 
be solved assuming the inverse of the submatrix in 
the upper right of the STM exists: 

A) = AAA \p f ~ ®ii(^o»^/)A)] 


Pf ~ ®2\(to’tf)Po T>22 (t( ) ,tf)p () 
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With these quantities and the STM, the primer 
vector at any time t can be solved for. 

Other important quantities include the 
Flamiltonian. While in general the expressions are 
different for the two-body system and the CRTBP 
system, they can be simplified to to the same 
expression, given below: 

H = -pv + pg (37) 

On an optimal arc this should be constant 0. 


KNITRO: NLP Code Usage 

All the solutions shown in this paper are found 
using the KNITRO nonlinear programming (NLP) 
algorithm. KNITRO, in its current form, is a 
package of three separate but not entirely distinct 
algorithms. They are, as named by KNITRO, an 
interior point conjugate gradient (CG) algorithm, an 
interior direct algorithm, and an SLP-EQP active set 
algorithm. 

Of these three, the interior direct method was the 
most successful and efficient in solving these 


problems. In fact, a rigorous study comparing 
KNITRO with another optimizer, SNOPT, resulted 
in an almost exclusive use of this algorithm. To 
explain the differences in general terms, the interior 
direct method is more or less an expansion of the 
interior CG that falls back into CG iterations as 
necessary. And essentially, the interior point CG 
algorithm “is a barrier method in which the 
subproblems are solved approximately by an SQP 
iteration with trust regions” (Byrd and Hribar 878). 

The main disadvantage of interior CG that 
motivates the design of the interior direct method is 
the use of a conjugate gradient iteration in the 
calculation of the tangential step. Basically, this 
iteration can be very expensive if the Hessian is ill- 
conditioned, in which case a direct linear algebra 
technique is better. Interior direct tackles this 
problem by introducing a line search method in 
place of this CG step. However, line searches are 
sometimes defective due to a non-convex quadratic 
model and rank deficiency of the Hessian. As a 
result, the interior direct method simply switches 
from a line search approach to the interior CG 
approach above as needed. This switch is 
advantageous because CG is uniquely equipped to 
deal with rank deficiencies due to its use of trust 
regions. 


Solution Cases 

The various solutions shown here will be 
presented progressively as the complexity in the type 
of problem being solved is increased. 

2-Body Orbit Transfer 

In order to verify the algorithm the first type of 
problem solved was a familiar orbit transfer in the 
two body frame. The following table shows the 
values for the two orbits involved in the transfer: 


Table 1: Two body Orbit Data 
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Here is the resulting solution: 


Orbital Trajectory: Converged Successfully [KNITRO Interior CG] 25 segments 



Figure 1: Two-body system orbit transfer solution 

The solution above was completed with 25 
segments. The initial guess is shown to be the black 
straight line between the dark blue and light blue 
orbits and the intermediate solutions that led to the 
final solution are shown in green. The final solution 
is red with each of the 25 individual segment nodes 
initial points indicated by squares. Here is a table for 
the final solution values: 


Table 2: Two body Orbit Transfer Results 


Parameter 

Symbol 

Value 

Units 

Departure 

t>o 

-67.065 

degrees 

Arrival 

u f 

181.782 

degrees 

Flight Time 

TOF 

29017.99 

seconds 

Initial Incr. 

Av 0 

.7828 

km/s 

Final Incr. 

A Vf 

1.744 

km/s 

Total Incr. 

A V t0 ,al 

2.527 

km/s 

The solution is 

locally 

optimal, as 

verified by 


Lawden’s necessary conditions. 

There are a number of things to note about this 
solution. First, the solution was obtained with the 
slightly modified approach using true anomalies 
instead of x’s. Second, while 25 segments are used 
and successfully converge, it does not appear that 25 
segments are the ideal number of segments to solve 
this problem. Fewer segments, if the initial 
conditions are tailored correctly is always more 
efficient. Third, while extensive knowledge of the 
problem does not seem to be necessary, some 
thought must go into the initial guess. In this case, 
the initial point on orbit 1 , the final point on orbit 2, 
and the transfer time are the independent variables 
that are used to create the straight line trajectory 
above. Choosing different points on the orbit seems 
to affect what local solution is converged to and how 
quickly this convergence happens, while changing 
the transfer time can even prevent convergence. 
Fourth, analytic gradients decrease computation time 
in this problem by more than 75% but the same 


number of major iterations are required. Finally, the 
optimizer spent most of its time optimizing after 
quickly finding a feasible solution. In this case only 
5 out of the 129 function evaluations were spent 
using a gravitational parameter other than the target 
value. 

CRTBP Moon Orbit to DRO Transfer 

Two impulse transfer trajectories are much more 
interesting in the circular restricted two body 
problem. The methodology described in this paper is 
general to any 3 -body system but in this instance the 
Earth-Moon system was chosen, partially because of 
the focus of the current Constellation project on 
returning to the moon. 


Table 3: CRTBP Moon Orbit and DRO Parameters 


Orbit 

Type 

x 0 , km 

y 0, rn/s 

T p , days 

Initial 

MO 

20000 

-495 

2.9375 

Final 

DRO 

100000 

633.62 

9.412 
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Figure 2: CRTBP DRO Transfer Solution 

The initial guess for the transfer is the series of 
solid black line segments in the above plot. Each 
segment is identical to its predecessor equidistant 
from the selected point on the initial orbit to the 
point on the final orbit. While the CRTBP is a much 
more complicated field in which to find impulsive 
trajectories, this solution was found relatively easily. 
The primer vector history shows that this an optimal 
two impulse transfer: 


i.i r 


Primer Vector vs Time: [KNITRO I nterior CG] 25 segments 



Figure 3: Primer Vector for DRO Solution 

It may not be apparent in the above plot the time 
rate of change of the primer vector is 0 at both the 
initial and final time. 

Table 4: CRTBP DRO Transfer Results 


Parameter 

Symbol 

Value 

Units 

Departure 

To 

7.02 

days 

Arrival 

Tf 

18.82 

days 

Flight Time 

TOF 

11.56 

days 

Initial Incr. 

Av 0 

273.26 

m/s 

Final Incr. 

Av f 

65.45 

m/s 

Total Incr. 

A V, ota i 

338.71 

m/s 


Again, like in the two-body orbit transfer, fewer 
segments can be used to find a convergence quicker, 
but comparing the work it takes to acquire these 
solutions reveals some drawbacks to using fewer 
segments. In the CRTBP system a zero gravitational 
parameter is no longer a force-free system. Thus 
each segment will not end up at the starting point of 
the next segment; in figure 2, each of the black 


segments along the solid line diverges noticeably 
from a straight line. The drawback is directly related 
to this: basically anytime an individual segment 
diverges significantly enough from the straight line 
by the dynamics the individual segments create 
wildly confusing guesses. By adding segments, the 
initial guess on transfer time can be worse, thus 
adding robustness. It is important to note that 
analytical gradients are a significant aid to reducing 
the solution time. 

CRTBP Earth Orbit to Lyapunov Orbit 

Transfers to libration point orbits such as 
Lyapunov and Halo orbits from initial Earth orbits 
can be difficult to design. Often optimal paths to 
these orbits are found by analyzing the unstable 
manifolds associated with them. These manifolds 
delineate paths that diverge quickly from the 
periodic orbits; if a spacecraft can hop on one in the 
vicinity of an Earth orbit the impulsive maneuver at 
the end of the trajectory is practically zero. The trick 
is finding a path that passes close to an initial 
parking orbit that requires a small impulse. 

If this algorithm is adept at picking up these 
manifold arcs and can reliably find transfers to these 
orbits without utilizing the dynamics associated with 
the manifold, it bodes well for the future 
development of a highly robust algorithm. If 
anything it could be useful for discovering new 
novel transfers. It appears that the algorithm may be 
sensitive enough, as the below converged transfer to 
a large Lyapunov orbit around Lj demonstrates. 
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Figure 4: Optimal Transfer to Lyapunov Orbit 


Table 5: Earth and Lyapunov Orbit Parameters 


Orbit 

Type 

x 0 , km 

To, km/s 

Initial 

EO 

75000 

2.305 

Final 

Lyapunov 

280000 

0.5902 


There are a number of qualities to note about 
this trajectory. First it appears to be a local 
minimum, as the primer vector plot helps confirm: 


Primer Vector vs Time: [KNITRO Interior CG] 25 segments 



Figure 5: Primer Vector for Lyapunov Transfer 

Although it is difficult to tell in figure 5, the rate 
is zero at the initial point. However, this is not true 
at the final point. It is speculated that this is due to 
the extremely small magnitude of the final 
maneuver. 

__Table_6^CRTBP^yapunov^ransferResults__ 

Parameter Symbol Value Units 


Departure 

to 

24.752 

days 

Arrival 

t f 

22.996 

days 

Flight Time 

TOF 

53.062 

days 

Initial Incr. 

Avo 

571.65 

m/s 

Final Incr. 

A Vf 

0.0045 

m/s 

Total Incr. 

A V to ,ai 

571.65 

m/s 


The initial guess was once again very important 
in finding a solution. The impact of the dynamics is 
much more obvious in the initial guess shown in 
figure 4. Especially near the Earth, the segments 
bend towards the Earth. Basically in order to get a 
solution the segments need to be long enough with a 
low enough initial velocity but not too long less they 
try to wrap around the primary. If this happens with 
the initial guess the segments split in a myriad of 
directions and only reconnect after many random 
passes around the system which is not likely to result 
in an optimal trajectory. 

An arbitrary set of orbits will not necessarily 
have a single pass optimal trajectory such as this. 
The two dimensional nature of this orbit transfer 
may also be a simpler trajectory to find. 

Perhaps the most interesting part about finding 
this solution is the need for a large number of 
segments. Even if the simple rules that have been 
established for having a good initial guess are 
followed as discussed in this paper, this trajectory 
was never found with less than 25 segments. Up 
until this point, fewer segments was always more 


advantageous, or seemingly so since run times could 
always be decreased. But it seems that having lots of 
segments not only makes the search more 
exploratory it may also keep the search well 
behaved. The evidence for this trend lies in the 
intermediate trajectories found from iteration to 
iteration. For fewer segments the intermediate 
trajectories can have very wild trajectory guesses. 
With 25 segments, the intermediate guesses usually 
stay disjointed longer but with a good initial guess 
the full trajectory rarely starts to make a huge jump 
between iterations to a radically different trajectory. 
The explanation might relate to the length of the 
segments. With only 1 long segment a small change 
in the initial state can make a large difference in the 
end point of that segment. But for each of several 
segments a small change has a much smaller effect 
on the entire trajectory. Up to this point the 
drawbacks for using many segments, mainly the 
increased run time, has outweighed the potential 
benefits But with analytical derivatives and a 
powerful NLP code, some of these drawbacks may 
be reduced significantly enough to make it work. 
This optimal trajectory is certainly a large step in 
that direction. 

Even with this solution, however, any transfer in 
the CRTBP frame appears to have less room for 
variation in the initial guess. This challenges the 
notion that it is truly robust with so much of the 
success of solving a problem depending on the 
character of the initial guess. 


CRTBP Earth Orbit to Halo Orbit 

While there seems to be a lot of success in 
finding transfers in the CRTBP frame to 2-d 
Lyapunov orbits, Halo orbits appear to present more 
trouble for the algorithm. Again, some of the 
transfers atempted may not have a good single pass 
transfer if the manifold does not pass closely to the 
orbit being transferred from. 
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Figure 6: Transfer to Halo Orbit, Top View 

However, while they may not be completely 
optimal, transfers to halo orbits that have very low 
final orbit transfer velocities can be found. 


Table 7: Earth Orbit and Halo Orbit Parameters 


Orbit 

Type 

x 0 , km 

z 0 , km 

V 0 , km/s 

Initial 

Final 

EO 

Halo 

90000 

374200 

0 

-27910 

2.305 

-.5344 


Trajectory in Rotating Frame: [KNITRO Interior CG] 25 segments 



Figure 7: Transfer to Halo Orbit, 3D View 


Table 8: Earth to Halo Orbit Transfer Results 


Parameter 

Symbol 

Value 

Units 

Departure 

To 

10.856 

days 

Arrival 

Tf 

10.994 

days 

Flight Time 

TOF 

19.711 

days 

Initial Incr. 

A v 0 

562.09 

m/s 

Final Incr. 

A Vf 

86.93 

m/s 

Total Incr. 

A V, ota i 

649.02 

m/s 


The primer vector plot reveals a lot about the 
optimality of the solution: 
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Figure 8: Primer Vector for Transfer to Halo 

Clearly the primer magnitude exceeds 1 and 
does so at three locations. This may be an optimal 
two-impulse solution that simply needs intermediate 
impulses. This presents an intriguing modification to 
the algorithm. If may be possible to modify the 
segment constraints so that the velocity increment 
between segments is simply added to the objective 
function rather than constrained to be continuous. 
While the analytic derivatives would need to be re- 
derived this would add capability to the algorithm to 
find multiple impulse solutions. 

CRTBP Halo Orbit to Halo Orbit 

Using primer vector theory to find orbit transfers 
is definitely not new. In fact, many of the concepts 
applied in this paper were used to find optimal 
trajectories in the CRTBP system in the available 
literature. 

In one, transfers between Halo orbits were 
found. These relatively simple CRTBP transfers 
were likewise modeled here to help verify 
implementation of the theory. Without a detailed 
account of the trajectories being analyzed it is hard 
to directly reproduce the results but it appears at 
least the character of the solution is similar: 


Trajectory inj^^ting Frame: [KNITRO Interior CG] 25 segments 



Figure 9: Halo to Halo Transfer 


Primer Vector vs Time: [KNITRO Interior CG] 25 segments 



Figure 10: Primer Vector for Halo to Halo Transfer 

The primer vector plot suggests that the initial 
and final times may need to shift earlier and later 
respectively as the primer rate does not go to zero at 
the end points. The fact that it converged may mean 
that the orbits are too sensitive. One remedy might 
be scaling. 


Conclusions 

The work presented here demonstrates that a 
new algorithm, a multiple direct shooting method 
that uses the gravitational parameter as an 
optimization variable, can be used to solve various 
orbital trajectory optimization problems. Some 
problems respond better than others to the algorithm. 
In the CRTBP system optimal transfers to periodic 
libration point orbits were the most surprising finds. 


Future modifications to the program may increase 
robustness further. This may be accomplished by 
including more advanced NLP codes, by adding the 
capability to handle intermediate impulses, or by 
designing an automated initial guess program. If the 
drawbacks can be ironed out, a robust algorithm may 
be on the horizon. 
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